function [ x, F, info]=Newton(fun, x0, opts, varargin)
%牛顿法  求min{F(x)}.x为n维向量, F为标量函数。
%F和F的梯度可由下面形式的matlab函数计算
%            function [F, g] = fun(x, p1, p2, ...)
% p1,p2,... 是这个函数的参数。
%
%  Call
%    [x, F, info] = Newton(fun, x0, opts)
%    [x, F, info] = Newton(fun, x0, opts,p1, p2, ...)
%
% 输入参数：
% fun  :  函数句柄。
% x0   :  初始点x0.
% opts :  最多五项的向量
%         opts(1)   :  Choice of method:  
%                      opts(1) = 0 : 阻尼牛顿法
%                      opts(1) = 1 : 混合（修正）牛顿法
%                      opts(1) = 2 : 拟牛顿法SR1
%                      opts(1) = 3 : 拟牛顿法BFGS
%                      otherwise   : 拟牛顿法DFP  (Default).
%         opts(2:4) :  终止准则参数。
%                      默认值: [1e-8  1e-8  1e-8]
%                      当不是混合牛顿法时只有opts(2)有效
%         opts(5)   :  最大迭代次数
%                      默认值: 10000.
%         opts(6)   :  单次线搜索的最大迭代次数
%                      默认值: 10.
%         如果输入参数不足6个，就取默认值。
% p1,p2,..  直接传给fun
%
% 输出参数：
% x     :  最优解
% F     :  最优值
% info  :  长度为3的向量
%          info(1) >  0 : 成功求解。
%                  =  0 : 达到最大迭代次数
%                  = -1 : x0 不是实向量
%                  = -2 : 迭代无法继续。
%                  = -3 : 线搜索出错。
%          info(2) = 函数调用次数
%          info(3) = 迭代次数
    [stop , n ] =check(x0);
    if stop, error('x must be a real valued vector'), end
    if  nargin < 3 | isempty(opts),  opts = 4; end
    if  opts(1) == 0,  
        opts = checkopts(opts, [0 1e-8 1e-8 1e-8 10000 10]);
    else,     
        opts = checkopts(opts, [4 1e-8 1e-8 1e-8 10000 10]);
    end
    method = opts(1);
    maxstep = opts(5);
    x = x0;F=fun(x,varargin{:});info = [1 0 0];
    lopts = [1  1e-3  1e-3  opts(6) 10];     % 线搜索参数。
                                             % lopts(1)=0精确，其他
                                             % 非精确
    % 阻尼newton法：
    if method == 0,
        Eps = opts(2);
        for i =1:maxstep,
            [F,g,G] = fun(x,varargin{:});
            info(2) = info(2)+1;
            if norm(g) <= Eps,
                info(3) = i;
                return
            end
            if cond(G) < 1000,
                d = -G\g;
                [x, F, g,infot]=linesearch(fun, x, F, g, d, lopts,varargin{:});
                info(2) = info(2)+infot(3);
                if infot(1)<=0
                    info(1) = -3;
                    info(3) = i;
                    %infot
                    return
                end
            else
                info(1) = -2;
                return
            end
        end
        if i == maxstep,
            info(1) = 0;
            info(3) = i;
        end
        return
    end
    % 混合newton法：
    if method == 1,
        Eps1 = opts(2);
        Eps2 = opts(3);
        Eps3 = opts(4);
        for i =1:maxstep,
            [F,g,G] = fun(x,varargin{:});
            info(2) = info(2)+1;
            if norm(g) <= Eps3,
                info(1) = 1;
                info(3) = i;
                return
            end
            if cond(G) < 1000,
                d = -G\g;
                dd = dot(d,g);
                nn = norm(d)*norm(g);
                if dd>Eps1*nn,
                    d = -d;
                else 
                    if abs(dd)<=Eps2*nn,
                        d = -g;
                    end
                end                
            else
                d = -g;
            end
            [x, F, g,infot]=linesearch(fun, x, F, g, d, lopts,varargin{:});
            info(2) = info(2)+infot(3);
            if infot(1)<=0
                [x, F, g,infot]=linesearch(fun, x, F, g, -g, lopts,varargin{:});
                info(2) = info(2)+infot(3);
                continue;
                info(1) = -3;
                info(3) = i;                
                %infot
                return
            end
        end
        if i == maxstep,
            info(1) = 0;
            info(3) = i;
        end
        return
    end
    %拟牛顿SR1：
    if method == 2,
        x = x0(:);Hp = eye(n);Hn=eye(n);Eps = opts(2);
        for i = 1:maxstep,
            [F,g] = fun(x,varargin{:});
            info(2) = info(2)+1;
            H = Hp;
            Hp = Hn;
            Hn = H;
            if norm(g) <= Eps,
                info(1) = 1;
                info(3) = i;
                return
            end
            d = -Hp*g;
            xp = x;
            gp = g;
            [x, F, g,infot]=linesearch(fun, xp, F, gp, d, lopts,varargin{:});
            info(2) = info(2)+infot(3);
            if infot(1)<=0
                Hp = eye(n); Hn = eye(n);
                continue;
                %info(1) = -3;
                %info(3) = i;
                %infot
                %return
            end
            s = x - xp;
            s = s(:);
            y = g - gp;
            y = y(:);
            shy = s - Hp*y;
            Hn = Hp + shy*shy'/dot(shy,y);
        end
        if i == maxstep,
            info(1) = 0;
            info(3) = i;
        end
    end    
    %拟牛顿BFGS：
    if method == 3,
        x = x0(:);Hp = eye(n);Hn=eye(n);Eps = opts(2);
        for i = 1:maxstep,
            [F,g] = fun(x,varargin{:});
            info(2) = info(2)+1;
            H = Hp;
            Hp = Hn;
            Hn = H;
            if norm(g) <= Eps,
                info(1) = 1;
                info(3) = i;
                return
            end
            d = -Hp*g;
            xp = x;
            gp = g;
            [x, F, g,infot]=linesearch(fun, xp, F, gp, d, lopts,varargin{:});
            info(2) = info(2)+infot(3);
            if infot(1)<=0
                Hp = eye(n); Hn = eye(n);
                continue;

                info(1) = -3;
                info(3) = i;
                %return
            end
            s = x - xp;
            s = s(:);
            y = g - gp;
            y = y(:);
            ys = dot(y,s);
            sy = s*y';
            Hn = Hp +(1+y'*Hp*y/ys)*(s*s'/ys)-(sy*Hp+Hp*sy')/ys;
        end
        if i == maxstep,
            info(1) = 0;
            info(3) = i;
        end
    end  
    %拟牛顿DFP：
    if method == 4,
        x = x0(:);Hp = eye(n);Hn=eye(n);Eps = opts(2);
        for i = 1:maxstep,
            [F,g] = fun(x,varargin{:});
            info(2) = info(2)+1;
            H = Hp;
            Hp = Hn;
            Hn = H;
            if norm(g) <= Eps,
                info(1) = 1;
                info(3) = i;
                return
            end
            d = -Hp*g;
            xp = x;
            gp = g;
            [x, F, g,infot]=linesearch(fun, xp, F, gp, d, lopts,varargin{:});
            info(2) = info(2)+infot(3);
            if infot(1)<=0
                Hp = eye(n); Hn = eye(n);
                continue;
                
                info(1) = -3;
                info(3) = i;
                %infot
                return
            end
            s = x - xp;
            s = s(:);
            y = g - gp;
            y = y(:);
            Hn = Hp + s*s'/dot(s,y) - Hp*y*y' *Hp/(y'*Hp*y);
        end
        if i == maxstep,
            info(1) = 0;
            info(3) = i;
        end
    end    
end
%============  Auxiliary functions  ========================

function [err,n] = check(x)
% CHECK - check x
%   
    err = 0; sx = size(x); n = max(sx);
    if  (min(sx) ~= 1) | ~isreal(x) | any(isnan(x(:))) | isinf(norm(x(:))) 
        err = -1; 
    end    
end
